library(tidyverse)
library(DESeq2)
library(glue)
library(DT)
library(BiocParallel)
library(patchwork)
library(writexl)
library(tibble)
library(clusterProfiler)
library(enrichplot)
library(org.Mm.eg.db)
library(AnnotationDbi)
library(ggpubr)
homedir <- "/central/groups/MazmanianLab/joeB"
wkdir <- glue("{homedir}/RNASEQ/fprau_bulk")
input_dir <- glue("{wkdir}/data/input")
interim_dir <- glue("{wkdir}/data/interim")
results_dir <- glue("{wkdir}/data/results")
fig_dir <- glue("{wkdir}/figures")
BiocParallel::register(BiocParallel::MulticoreParam(6))ASO F. Prausnitzii Treatment Bulk RNA-Seq Analysis
Overview
- Imports raw gene expression data and metadata into DESeq2 objects
- Performs differential expression analysis
- Visualizes results
- Performs enrichment analysis
- Saves results tables
Methods
See exact packages below for versions. At a high level we use DESeq2 for differential expression analysis, clusterProfiler for enrichment analysis, and ggplot2 for visualization.
Analysis Setup
Environment setup
Prepping DESeq2 object
# Loading metadata
meta <- read_csv(glue("{input_dir}/metadata_BulkRNASeq_Samples_AM.csv")) %>%
janitor::clean_names() %>%
mutate(group = glue("{genotype} {treatment}")) %>%
mutate(group = factor(group, levels = c("ASO Vehicle", "WT Vehicle", "ASO Fp"))) %>%
mutate_at(vars(extraction_date, collection_date), lubridate::mdy) %>%
glimpse()
# Examing the metadata
DT::datatable(meta)MID BRAIN
- extraction date all the same
- cage coincides w/ conditions, - not a helpful covar in this case
- collection date is even across conditions - useful var to check for batch effects
- all mice are from the same cohort
LARGE INTESTINE
- cohort, cage, extraction date, collection date all might be useful covars
# Loading gene expression data
gex <- read_csv(glue("{input_dir}/20250416_ParkinsonsAVW_rawcounts.csv")) %>%
glimpse()
gene_map <- gex %>% dplyr::select(gene_id, gene_name, gene_type) %>% distinct()
gene_mat <- gex %>%
column_to_rownames(var = "gene_id") %>%
select(-c(gene_name, gene_type)) %>%
as.matrix()
meta_gi <- meta %>% filter(tissue == "Large intestine") %>%
mutate_at(vars(genotype, treatment, cohort, cage_number), as.factor)
rownames(meta_gi) <- meta_gi$sample_id
meta_mb <- meta %>% filter(tissue == "Mid Brain") %>%
mutate_at(vars(genotype, treatment, cohort, cage_number), as.factor)
rownames(meta_mb) <- meta_mb$sample_id
# Creating DESeq2 object
# DeSeq2 is unable to fit the glm when more than one factor is included in the design
dds_mb <- DESeqDataSetFromMatrix(countData = gene_mat[,rownames(meta_mb)],
colData = meta_mb,
design = ~ group
)
dds_gi <- DESeqDataSetFromMatrix(countData = gene_mat[,rownames(meta_gi)],
colData = meta_gi,
design = ~ group
)Minor pre-filtering - removing genes with total sum across samples <= 10 counts
nrow(dds_mb)
dds_mb <- dds_mb[rowSums(counts(dds_mb)) >= 10, ]
nrow(dds_mb)
# 78298 to 47018 genes remaining
nrow(dds_gi)
dds_gi <- dds_gi[rowSums(counts(dds_gi)) >= 10, ]
nrow(dds_gi)
# 78298 to 41317 genes remainingDifferential Expression Analysis
Running DESeq2 will perform the following key steps:
Estimation of size factors
Estimation of dispersion
Negative Binomial GLM fitting and Wald statistic
dds_mb <- DESeq(dds_mb)
dds_mb
dds_gi <- DESeq(dds_gi)
dds_gicomps <- tribble(
~numerator, ~denominator,
"WT Vehicle", "ASO Vehicle",
"ASO Fp", "WT Vehicle",
"ASO Fp", "ASO Vehicle"
)
res_mb <- list()
res_gi <- list()
for (i in 1:nrow(comps)) {
vs <- glue("{comps$numerator[i]}_vs_{comps$denominator[i]}") %>% gsub(" ", "_", .)
message(glue(
"Calculating differential expression for ",
"{comps$numerator[i]} vs {comps$denominator[i]}"
))
res_mb[[vs]] <- results(
dds_mb,
contrast = c("group", comps$numerator[i], comps$denominator[i])
)
res_gi[[vs]] <- results(
dds_gi,
contrast = c("group", comps$numerator[i], comps$denominator[i])
)
}
saveRDS(dds_mb, glue("{interim_dir}/deseq2_mb_{Sys.Date()}.rds"))
saveRDS(dds_gi, glue("{interim_dir}/deseq2_gi_{Sys.Date()}.rds"))
saveRDS(res_mb, glue("{interim_dir}/deseq2_mb_res_{Sys.Date()}.rds"))
saveRDS(res_gi, glue("{interim_dir}/deseq2_gi_res_{Sys.Date()}.rds"))Visualizing results
Reading in processed results
dds_mb <- readRDS(glue("{interim_dir}/deseq2_mb_2025-05-09.rds"))
dds_gi <- readRDS(glue("{interim_dir}/deseq2_gi_2025-05-09.rds"))
res_mb <- readRDS(glue("{interim_dir}/deseq2_mb_res_2025-05-09.rds"))
res_gi <- readRDS(glue("{interim_dir}/deseq2_gi_res_2025-05-09.rds"))Visualzing p-value histograms
plotting funcs
plot_pval_hist <- function(res) {
title = res@elementMetadata$description[5]
res %>%
as.data.frame() %>%
ggplot(aes(x = pvalue)) +
geom_histogram(binwidth = 0.01) +
theme_minimal() +
ggtitle(title)
}
plot_volcano <- function(res) {
title <- res@elementMetadata$description[5]
title <- res@elementMetadata$description[5] %>%
strex::str_after_first("group ")
numerator_group <- title %>% strex::str_before_first(" vs ")
denominator_group <- title %>% strex::str_after_first(" vs ")
color_pal <- c("#EF433D", "#06A0DD", "gray80")
names(color_pal) <- c(
glue("{numerator_group} Up"),
glue("{denominator_group} Up"),
"ns"
)
# Convert to dataframe and add gene symbols
plot_df <- res %>%
as.data.frame() %>%
drop_na(padj) %>%
rownames_to_column(var = "gene_id") %>%
left_join(gene_map, by = "gene_id") %>%
mutate(sig = case_when(
padj <= 0.05 & log2FoldChange > 0.2 ~ glue("{numerator_group} Up"),
padj <= 0.05 & log2FoldChange < -0.2 ~ glue("{denominator_group} Up"),
TRUE ~ "ns"
))
# Identify top 10 genes by significance (lowest padj)
top_genes <- plot_df %>%
filter(sig !="ns") %>%
filter(padj < 0.05) %>%
arrange(padj) %>%
head(15)
# Create the plot with labels
ggplot(plot_df, aes(x = log2FoldChange, y = -log10(padj), color = sig)) +
geom_point(alpha = 0.75) +
theme_bw() +
ggtitle(title) +
scale_color_manual(values = color_pal) +
theme(legend.position = "bottom") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(color = NULL, y = expression(-log[10](padj)), x = expression(log[2](Fold~Change))) +
# Add ggrepel labels for top genes
ggrepel::geom_text_repel(
data = top_genes,
aes(label = gene_name),
size = 3,
max.overlaps = Inf,
force = 50,
force_pull = 0,
nudge_x = ifelse(top_genes$log2FoldChange > 0, 1, -1),
nudge_y = 1,
box.padding = 0.8,
point.padding = 0.5,
segment.color = "grey",
segment.size = 0.2,
segment.alpha = 0.75,
direction = "both",
# Add these parameters for curved segments:
segment.curvature = 0.3, # Controls how curved the lines are (0 = straight)
segment.ncp = 3, # Number of control points (higher = smoother curves)
segment.angle = 10 # Angle of the curve
)
}P-value Distributions
p_hist_mb <- res_mb %>%
purrr::map(~plot_pval_hist(.x))
p_hist_gi <- res_gi %>%
purrr::map(~plot_pval_hist(.x))
p_hist_mb_patch <- p_hist_mb %>% purrr::reduce(`/`) + plot_layout(ncol = 1)
p_hist_gi_patch <- p_hist_gi %>% purrr::reduce(`/`) + plot_layout(ncol = 1)
ggsave(p_hist_mb_patch,
filename = glue("{fig_dir}/deseq2_mb_pval_hist.png"), width = 6, height = 10)
ggsave(p_hist_gi_patch,
filename = glue("{fig_dir}/deseq2_gi_pval_hist.png"), width = 6, height = 10)Takeaways - P-value Distributions
- P-value distributions are robust for Large Intestine, but not for Mid Brain - which seem to point towards a null hypothesis.
Dimension Reduction - PCA
vsd_mb <- vst(dds_mb, blind=FALSE)
vsd_gi <- vst(dds_gi, blind=FALSE)
p_pca_mb <- plotPCA(vsd_mb, intgroup=c("group")) +
theme_bw() +
coord_fixed() +
ggtitle("Mid Brain") +
scale_color_manual(
values = c("ASO Vehicle" = "#A6CEE3", "ASO Fp" = "#1F78B4", "WT Vehicle" = "#B2DF8A"),
limits = c("ASO Vehicle", "ASO Fp", "WT Vehicle")
)
p_pca_gi <- plotPCA(vsd_gi, intgroup=c("group")) +
theme_bw() +
coord_fixed() +
ggtitle("Large Intestine") +
scale_color_manual(
values = c("ASO Vehicle" = "#A6CEE3", "ASO Fp" = "#1F78B4", "WT Vehicle" = "#B2DF8A"),
limits = c("ASO Vehicle", "ASO Fp", "WT Vehicle")
)
ggsave(p_pca_mb,
filename = glue("{fig_dir}/deseq2_mb_pca.png"), width = 6, height = 5)
ggsave(p_pca_gi,
filename = glue("{fig_dir}/deseq2_gi_pca.png"), width = 6, height = 5)
Takeaways - PCA
- PCA plots show some separation between conditions, but differences are modest. There is substantial similarity between a handful of some samples across conditions in each tissue.
Volcano Plots
volcano_figs <- glue("{fig_dir}/DESeq2_VolcanoPlots")
dir.create(volcano_figs, showWarnings = FALSE, recursive = TRUE)
for (comp in names(res_gi)) {
p_vol_gi <- res_gi[[comp]] %>% plot_volcano()
p_vol_mb <- res_mb[[comp]] %>% plot_volcano()
ggsave(
p_vol_gi,
filename = glue("{volcano_figs}/GI_volcano_{comp}.png"),
width = 5,
height = 5
)
ggsave(
p_vol_mb,
filename = glue("{volcano_figs}/MB_volcano_{comp}.png"),
width = 5,
height = 5
)
}Visualzing similarities in gene expression shifts between conditions
names(res_gi)
wt_v_asoc <- res_gi[["WT_Vehicle_vs_ASO_Vehicle"]] %>%
as.data.frame() %>%
mutate(sig = case_when(
padj <= 0.05 & log2FoldChange > 0.2 ~ "Up",
padj <= 0.05 & log2FoldChange < -0.2 ~ "Down",
TRUE ~ "ns"
)) %>%
dplyr::rename_all( ~ paste0(., "__wt_v_asoc")) %>%
rownames_to_column(var = "gene_id")
asofp_v_asoc <- res_gi[["ASO_Fp_vs_ASO_Vehicle"]] %>%
as.data.frame() %>%
mutate(sig = case_when(
padj <= 0.05 & log2FoldChange > 0.2 ~ "Up",
padj <= 0.05 & log2FoldChange < -0.2 ~ "Down",
TRUE ~ "ns"
)) %>%
dplyr::rename_all( ~ paste0(., "__asofp_v_asoc")) %>%
rownames_to_column(var = "gene_id")
df_comb <- full_join(wt_v_asoc, asofp_v_asoc) %>%
left_join(gene_map, by = "gene_id") %>%
mutate(sig_comb = case_when(
sig__wt_v_asoc == "Up" & sig__asofp_v_asoc == "Up" ~ "WT & ASO Fp Up",
sig__wt_v_asoc == "Down" & sig__asofp_v_asoc == "Down" ~ "WT & ASO Fp Down",
TRUE ~ "Mixed"
)) %>%
drop_na(pvalue__wt_v_asoc, pvalue__asofp_v_asoc) %>%
glimpse()
df_comb$sig_comb %>% table()
top_genes <- df_comb %>%
filter(sig_comb != "Mixed")
df_comb %>% glimpse()
color_pal <- c("#EF433D", "#06A0DD", "gray80")
names(color_pal) <- c("WT & ASO Fp Up", "WT & ASO Fp Down", "Mixed")
set.seed(123)
p_comb <- df_comb %>%
drop_na(log2FoldChange__wt_v_asoc, log2FoldChange__asofp_v_asoc) %>%
# Reorder the factor levels to make "Mixed" appear first (below other points)
mutate(sig_comb = factor(sig_comb, levels = c("Mixed", "WT & ASO Fp Up", "WT & ASO Fp Down"))) %>%
arrange(sig_comb) %>%
ggplot(aes(x = log2FoldChange__wt_v_asoc, y = log2FoldChange__asofp_v_asoc)) +
geom_point(aes(color = sig_comb), alpha = 0.75) +
# geom_abline(slope = 1, linetype = "dotted") + # y ~ x line
# geom_abline(slope = -1, linetype = "dotted") + # y ~ -x line
scale_color_manual(values = color_pal) +
labs(
color = NULL,
x = expression(log[2](Fold~Change[WT~Vehicle~vs~ASO~Vehicle])),
y = expression(log[2](Fold~Change[ASO~Fp~vs~ASO~Vehicle])),
color = NULL
) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
ggpubr::stat_cor(
aes(label = paste(..r.label.., ..p.label.., sep = "~`,`~")),
method = "pearson", label.x = -5, label.y = 5, na.rm = TRUE
) +
ggrepel::geom_label_repel(
data = top_genes,
aes(label = gene_name, color = sig_comb),
size = 3,
max.overlaps = Inf,
force = 5,
force_pull = 0,
nudge_x = ifelse(top_genes$log2FoldChange > 0, 1, -1),
nudge_y = 1,
box.padding = 0.8,
point.padding = 0.5,
segment.color = "black",
segment.size = 0.2,
segment.alpha = 0.6,
direction = "both",
# Add these parameters for curved segments:
segment.curvature = 0.3, # Controls how curved the lines are (0 = straight)
segment.ncp = 3, # Number of control points (higher = smoother curves)
segment.angle = 10 # Angle of the curve
) +
coord_fixed() +
theme_bw()
ggsave(
p_comb,
filename = glue("{fig_dir}/cross_condition_L2fc_scatter.png"),
width = 7,
height = 5
)Takeaways - Cross Condition L2FC Scatter
- Cross condition L2FC scatter plot highlights a handful of genes that are enriched in two comparisons of interest in the same direction. 1 - Up in WT compared to ASO Vehicle, 2 - Up in ASO Fp compared to ASO Vehicle. We correlated the log2fc of all genes retained in the DESeq2 analysis and observe a significant positive correlation between the two conditions - suggesting that F. prausnitzii treatment in ASO shifts the large intestine expression profile towards the WT Vehicle profile with respect to an altered ASO Vehicle background. Points that are colored pass significance in both comparisons are highlighted with a label.
Enrichment Analysis
source(glue("{wkdir}/notebook/R_scripts/quick_enrich.R"))
gofigs <- glue("{fig_dir}/GO_enrichment")
godata <- glue("{interim_dir}/GO_enrichment")
# dir.create(godata, showWarnings = FALSE, recursive = TRUE)
for (ont_type in c("BP")) {
for (comp in names(res_gi)) {
message(glue("Running enrichment analysis for {comp}"))
enrich_out_gi <- run_GO_enrich(
res = res_gi[[comp]],
alpha = 0.1,
lfc = 0.2,
ont = ont_type,
top_n = 20,
plot_file = glue("{gofigs}/GO_dotplot_GI_{comp}_{ont_type}.png")
)
saveRDS(enrich_out_gi, glue("{godata}/enrich_out_GI_{comp}_{ont_type}_{Sys.Date()}.rds"))
enrich_out_mb <- run_GO_enrich(
res = res_mb[[comp]],
alpha = 0.1,
lfc = 0.2,
ont = ont_type,
top_n = 20,
plot_file = glue("{gofigs}/GO_dotplot_MB_{comp}_{ont_type}.png")
)
saveRDS(enrich_out_mb, glue("{godata}/enrich_out_MB_{comp}_{ont_type}_{Sys.Date()}.rds"))
}
}Large Intestine
Mid Brain
Takeaways - GO Enrichment
GO enrichment analysis shows robust hits for the Large Intestine, but not for the Mid Brain.
GO analysis of the large intenstine show some interesting results -
- the ASO Fp vs ASO Vehicle contrast (left panel) is enriched for actin-filament organisation, wound healing, skeletal morphogenesis and cartilage development, pointing to cytoskeletal re-arrangement and tissue-repair pathways rather than immune activation. This suggests that the F. prausnitzii ASO is inducing a “healing” response.
- Immune pathways dominate whenever the WT background is involved – both contrasts that include WT animals (centre and right panels) rank lymphocyte activation/differentiation, leukocyte migration, and regulation of immune-system processes at the very top, indicating that genotype (WT vs ASO) is primarily distinguished by immune gene programmes.
- The WT Vehicle vs ASO Vehicle comparison shows the strongest proportional shift – its GeneRatio values climb to ~0.09 (higher than the ≤ 0.06 seen in the other two plots), suggesting that baseline genetic difference out-weighs Fp treatment in terms of the fraction of genes hitting each GO term.
- Cell–cell adhesion appears as a shared immune theme, absent from the Fp-only contrast – terms such as “regulation of cell-cell adhesion” and “leukocyte cell–cell adhesion” recur in the two immune-heavy panels but not in the Fp-vehicle panel, underscoring that adhesion/trafficking of immune cells is a WT/ASO signature rather than a direct consequence of Fp treatment.
Saving Results Tables
# Convert DESeq2 results to data frame
format_results_df <- function(res) {
df <- res %>%
as.data.frame() %>%
rownames_to_column(var = "gene_id") %>%
left_join(gene_map, by = "gene_id") %>%
glimpse()
return(df)
}
gi_res_list <- res_gi %>% purrr::map(~format_results_df(.x))
mb_res_list <- res_mb %>% purrr::map(~format_results_df(.x))
# Write to Excel file
write_xlsx(gi_res_list, path = glue("{results_dir}/DESeq2_GI_results.xlsx"))
write_xlsx(mb_res_list, path = glue("{results_dir}/DESeq2_MB_results.xlsx"))
# Saving GO enrichment results
godata_files <- list.files(godata, full.names = TRUE, pattern = "2025-05-10.rds")
godata_list <- godata_files %>%
purrr::set_names(basename(.) %>% gsub("enrich_out_|_2025-05-10.rds", "", .)) %>%
purrr::map(~readRDS(.x)$results %>% as.data.frame())
# Convert list of lists to a single data frame
write_xlsx(godata_list, path = glue("{results_dir}/GO_enrichment_results.xlsx"))Takeaways - ALL
- P-value distributions are robust for Large Intestine, but not for Mid Brain - which seem to point towards a null hypothesis.
- PCA plots show some separation between conditions, but differences are modest. There is substantial similarity between a handful of some samples across conditions in each tissue.
- Volcano plots show robust differences between groups of interest in primarily the Large Intestine. In the Mid Brain, there are fewer differences, and they seem primarily related to the genetic construct of the ASO/WT background.
- Cross condition L2FC scatter plot highlights a handful of genes that are enriched in two comparisons of interest in the same direction. 1 - Up in WT compared to ASO Vehicle, 2 - Up in ASO Fp compared to ASO Vehicle. We correlated the log2fc of all genes retained in the DESeq2 analysis and observe a significant positive correlation between the two conditions - suggesting that F. prausnitzii treatment in ASO shifts the large intestine expression profile towards the WT Vehicle profile with respect to an altered ASO Vehicle background.
- GO enrichment analysis shows robust hits for the Large Intestine, but not for the Mid Brain.
- GO analysis of the large intenstine show some interesting results -
- the ASO Fp vs ASO Vehicle contrast (left panel) is enriched for actin-filament organisation, wound healing, skeletal morphogenesis and cartilage development, pointing to cytoskeletal re-arrangement and tissue-repair pathways rather than immune activation. This suggests that the F. prausnitzii ASO is inducing a “healing” response.
- Immune pathways dominate whenever the WT background is involved – both contrasts that include WT animals (centre and right panels) rank lymphocyte activation/differentiation, leukocyte migration, and regulation of immune-system processes at the very top, indicating that genotype (WT vs ASO) is primarily distinguished by immune gene programmes.
- The WT Vehicle vs ASO Vehicle comparison shows the strongest proportional shift – its GeneRatio values climb to ~0.09 (higher than the ≤ 0.06 seen in the other two plots), suggesting that baseline genetic difference out-weighs Fp treatment in terms of the fraction of genes hitting each GO term.
- Cell–cell adhesion appears as a shared immune theme, absent from the Fp-only contrast – terms such as “regulation of cell-cell adhesion” and “leukocyte cell–cell adhesion” recur in the two immune-heavy panels but not in the Fp-vehicle panel, underscoring that adhesion/trafficking of immune cells is a WT/ASO signature rather than a direct consequence of Fp treatment.
Environment Info
sessionInfo()R version 4.2.0 (2022-04-22)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 9.3 (Plow)
Matrix products: default
BLAS/LAPACK: /central/groups/MazmanianLab/joeB/software/mambaforge/envs/pdmbsR/lib/libopenblasp-r0.3.23.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] ggpubr_0.6.0 org.Mm.eg.db_3.16.0
[3] AnnotationDbi_1.60.2 enrichplot_1.18.0
[5] clusterProfiler_4.6.0 writexl_1.5.0
[7] patchwork_1.1.2 BiocParallel_1.32.6
[9] DT_0.28 glue_1.6.2
[11] DESeq2_1.38.3 SummarizedExperiment_1.28.0
[13] Biobase_2.58.0 MatrixGenerics_1.10.0
[15] matrixStats_1.0.0 GenomicRanges_1.50.2
[17] GenomeInfoDb_1.34.9 IRanges_2.32.0
[19] S4Vectors_0.36.2 BiocGenerics_0.44.0
[21] lubridate_1.9.2 forcats_1.0.0
[23] stringr_1.5.0 dplyr_1.1.2
[25] purrr_1.0.1 readr_2.1.4
[27] tidyr_1.3.0 tibble_3.2.1
[29] ggplot2_3.4.4 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] shadowtext_0.1.3 backports_1.4.1 fastmatch_1.1-3
[4] plyr_1.8.8 igraph_1.5.0 lazyeval_0.2.2
[7] splines_4.2.0 digest_0.6.30 yulab.utils_0.0.6
[10] htmltools_0.5.5 GOSemSim_2.24.0 viridis_0.6.3
[13] GO.db_3.16.0 fansi_1.0.3 magrittr_2.0.3
[16] memoise_2.0.1 tzdb_0.4.0 Biostrings_2.66.0
[19] annotate_1.76.0 graphlayouts_0.8.0 timechange_0.2.0
[22] colorspace_2.1-0 blob_1.2.4 ggrepel_0.9.3
[25] xfun_0.39 crayon_1.5.2 RCurl_1.98-1.12
[28] jsonlite_1.8.3 scatterpie_0.2.3 ape_5.7-1
[31] polyclip_1.10-4 gtable_0.3.3 zlibbioc_1.44.0
[34] XVector_0.38.0 DelayedArray_0.24.0 car_3.1-2
[37] abind_1.4-5 scales_1.2.1 DOSE_3.24.0
[40] DBI_1.1.3 rstatix_0.7.2 Rcpp_1.0.10
[43] viridisLite_0.4.2 xtable_1.8-4 gridGraphics_0.5-1
[46] tidytree_0.4.2 bit_4.0.5 htmlwidgets_1.6.2
[49] httr_1.4.6 fgsea_1.24.0 RColorBrewer_1.1-3
[52] pkgconfig_2.0.3 XML_3.99-0.9 farver_2.1.1
[55] locfit_1.5-9.8 utf8_1.2.2 ggplotify_0.1.0
[58] tidyselect_1.2.0 rlang_1.1.1 reshape2_1.4.4
[61] munsell_0.5.0 tools_4.2.0 cachem_1.0.8
[64] downloader_0.4 cli_3.6.1 generics_0.1.3
[67] RSQLite_2.3.1 gson_0.1.0 broom_1.0.5
[70] evaluate_0.21 fastmap_1.1.1 yaml_2.3.6
[73] ggtree_3.6.2 knitr_1.40 bit64_4.0.5
[76] tidygraph_1.3.0 KEGGREST_1.38.0 ggraph_2.1.0
[79] nlme_3.1-162 aplot_0.1.10 compiler_4.2.0
[82] png_0.1-8 ggsignif_0.6.4 treeio_1.22.0
[85] tweenr_2.0.2 geneplotter_1.76.0 stringi_1.7.8
[88] lattice_0.21-8 Matrix_1.6-3 vctrs_0.6.3
[91] pillar_1.9.0 lifecycle_1.0.3 data.table_1.14.8
[94] cowplot_1.1.1 bitops_1.0-7 qvalue_2.30.0
[97] R6_2.5.1 gridExtra_2.3 codetools_0.2-19
[100] MASS_7.3-60 withr_2.5.0 GenomeInfoDbData_1.2.9
[103] parallel_4.2.0 hms_1.1.3 grid_4.2.0
[106] ggfun_0.0.9 HDO.db_0.99.1 rmarkdown_2.22
[109] carData_3.0-5 ggforce_0.4.1